function bvp

%  Solves the BVP:
%       y'' + p(x)y' + q(x)y= f(x)   for xL < x < xR  
%  where 
%      y(xL) = yL  and y(xR) = yR

clear *
clf

n=100;
ep=1;

hold on
% get(gcf)
set(gcf,'Position', [838 464 530 280]);

% set boundary conditions
xL=0; yL=0;
xR=1; yR=0;

% generate the points along the x-axis, x(1)=xL and x(n)=xR
x=linspace(xL,xR,n);
h=x(2)-x(1);

% calculate the coefficients of finite difference equation
a=zeros(1,n-2); b=zeros(1,n-2); c=zeros(1,n-2);
for i=1:n-2
	a(i) = -2+h*h*q(x(i+1),ep);
	b(i) = 1-0.5*h*p(x(i+1),ep);
	c(i) = 1+0.5*h*p(x(i+1),ep);
	f(i) = h*h*rhs(x(i+1));
end;
f(1) = f(1)-yL*b(1);
f(n-2) = f(n-2)-yR*c(n-2);
% solve the tri-diagonal matrix problem
y=tri(a,b,c,f);
y=[yL, y, yR];

plot(x,y,'-','LineWidth',1,'MarkerSize',7)

xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')

%axis([0 1 -0.04 0.04]);
%legend(' Exact',' N = 2',' N = 6',3)
box on
%set(gca,'ytick',[-0.04  -0.02 0 0.02 0.04]);
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
%set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
hold off

function g=rhs(x,ep)
g=-sin(2*pi*x);

function g=q(x,ep)
g=-1;

function g=p(x,ep)
g=0;

% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end
